import pandas as pd
import numpy as np
import os
import shap
import xgboost as xgb
from sklearn.linear_model import Ridge
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import StandardScaler
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import precision_score, recall_score, accuracy_score, roc_curve, \
roc_auc_score, precision_recall_curve, confusion_matrix, r2_score, mean_absolute_error, \
mean_squared_error
import pandas_profiling as pp
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:60% !important; }</style>"))
Read the data
!dir data
data_dir = 'data'
df = pd.read_csv(os.path.join(data_dir, 'MEPS_data_preprocessed.csv'))
df.reset_index(drop=True, inplace=True)
df_raw = df.copy()
df.head()
Change categorical variables to strings (even though they're represented as numbers)
cats = ['REGION','MARRY31X','EDRECODE','FTSTU31X','ACTDTY31','HONRDC31',
'RTHLTH31','MNHLTH31','HIBPDX','CHDDX','ANGIDX','MIDX','OHRTDX','STRKDX',
'EMPHDX','CHBRON31','CHOLDX','CANCERDX','DIABDX','JTPAIN31','ARTHDX',
'ARTHTYPE','ASTHDX','ADHDADDX','PREGNT31','WLKLIM31','ACTLIM31','SOCLIM31',
'COGLIM31','DFHEAR42','DFSEE42','ADSMOK42','PHQ242','EMPST31','POVCAT15','INSCOV15']
df = df.copy()
for col in cats:
df[col] = df[col].astype(str)
# # Show summary report
# df.profile_report(title='Summary Report')
Let's have a closer look at the predicted variable (HEALTHEXP)
df['HEALTHEXP'].hist(bins=100)
df['HEALTHEXP'].describe(percentiles=np.linspace(0,1,21))
health_exp = df['HEALTHEXP'].values
plt.plot(np.linspace(0,100,21), np.percentile(health_exp, np.linspace(0,100,21)),'o');
plt.axhline(y=np.mean(health_exp), color='red', label = 'mean')
plt.xlabel('Percentile');
plt.ylabel('Healthexp');
plt.legend();
plt.plot(np.linspace(0,100,21), np.percentile(health_exp, np.linspace(0,100,21)),'o');
plt.axhline(y=np.mean(health_exp), color='red', label = 'mean')
plt.xlabel('Percentile');
plt.ylabel('Healthexp');
plt.yscale('log')
plt.legend();
plt.figure(figsize=(6,6))
df[['HEALTHEXP']].boxplot();
plt.figure(figsize=(6,6))
plt.boxplot(df['HEALTHEXP']);
plt.title('Health Exp')
plt.yscale('log')
Notes on predicted variable:
- 14.3 % of observations are equal to 0.
- The distribution is positevely skewed with many high-value outliers (very high treatment costs).
- About 80 % of observations are below the mean. Std is three times higher than the mean.
cut_off = 100_000
df = df[df['HEALTHEXP']<cut_off]
Prepare features
# Drop panel number (not meant to be predictive) and sample weights
df.drop(columns = ['PANEL', 'PERSONWT'], inplace=True)
df.head()
y = df.pop('HEALTHEXP')
One-hot encoding for all the categorical features
cat_cols_one_hot = []
for col in cats:
var_one_hot = pd.get_dummies(df[col],prefix=col, drop_first = True)
df = df.drop(columns=[col])
df = pd.concat([df, var_one_hot], axis=1)
for cl in var_one_hot.columns:
cat_cols_one_hot.append(cl)
df.head()
More feature engineering in next HW ...
def train_print_save_models_reg(mods, dfTrain, dfTest, yTrain, yTest):
trained_models = {}
for model_name, mod_object in mods.items():
arTrain = np.asarray(dfTrain)
arTest = np.asarray(dfTest)
md = mod_object.fit(arTrain,yTrain)
y_pred_test = md.predict(arTest)
y_pred_test[y_pred_test<0]=0
y_pred_train = md.predict(arTrain)
mae_train = mean_absolute_error(yTrain, y_pred_train)
mae_test = mean_absolute_error(yTest, y_pred_test)
rmse_train = np.sqrt(mean_squared_error(yTrain, y_pred_train))
rmse_test = np.sqrt(mean_squared_error(yTest, y_pred_test))
r2_train = r2_score(yTrain, y_pred_train)
r2_test = r2_score(yTest, y_pred_test)
print(model_name)
print('Training R^2: {}, MAE:{}, RMSE:{} '.format(r2_train, mae_train, rmse_train))
print('Test R^2:{}, MAE:{}, RMSE:{}'.format(r2_test, mae_test, rmse_test) )
trained_models[model_name] = md
return trained_models
mods = {
'XGB':xgb.XGBRegressor(max_depth=3),
'LinReg':LinearRegression(),
'Ridge_Reg':Ridge(alpha=.5)
}
dfTrain, dfTest, yTrain, yTest = train_test_split(df, y, random_state=42)
len(dfTrain), len(dfTest)
models = train_print_save_models_reg(mods, dfTrain, dfTest, yTrain, yTest)
import tensorflow as tf
from keras.models import Sequential
from keras.optimizers import RMSprop
from keras.layers import Dense, Dropout, Embedding, Reshape
from keras.backend import clear_session
from keras.optimizers import Adam
from keras import backend as K
from keras.models import Model
from keras.layers import Dense
clear_session()
def r2(y_true, y_pred):
SS_res = K.sum(K.square( y_true-y_pred ))
SS_tot = K.sum(K.square( y_true - K.mean(y_true) ) )
return ( 1 - SS_res/(SS_tot + K.epsilon()) )
model = Sequential([
Dense(64, activation=tf.nn.relu, input_shape=[dfTrain.shape[1]]),
Dropout(0.3),
Dense(128, activation='relu'),
Dense(64, activation='relu'),
Dense(1)
])
optimizer = Adam(1e-4)
model.compile(loss='mae',
optimizer = optimizer,
metrics=[r2])
model.fit(x=dfTrain,
y=yTrain, epochs=1000,
batch_size =1024,
validation_data=(dfTest,yTest),
verbose=2,
validation_freq=1, use_multiprocessing=True)
Simple ANN
y_pred = model.predict(dfTest)
r2_score(yTest, y_pred)
mae = mean_absolute_error(yTest, y_pred)
rmse = np.sqrt(mean_squared_error(yTest, y_pred))
print('MAE: {}, RMSE: {}'.format(mae, rmse))
models['ANN'] = model
For some selected observation from this dataset, calculate the model predictions for model (1)¶
arTrain = np.asarray(dfTrain)
arTest = np.asarray(dfTest)
model1 = models['XGB']
model1
obs_num_1 = 10
arobs = np.reshape(arTest[obs_num_1], (1,-1))
dfTest.iloc[obs_num_1,:]
#wartość rzeczywista
yTest.iloc[obs_num_1]
#predykcja
model1.predict(arobs)[0]
For an observation selected in (2), calculate the decomposition of model prediction using LIME / live / lime / localModel or similar technique (packages for R: live, live, localModel, iml, packages for python: lime).¶
import lime
from lime.lime_tabular import LimeTabularExplainer as LTE
cat_features = np.asarray([i for (i,col) in enumerate(dfTrain.columns) if col in cat_cols_one_hot])
feat_names = np.asarray(list(dfTrain.columns))
dfTrain.columns[cat_features]
explainer = LTE(arTrain, mode='regression', feature_names=feat_names,
class_names=['asd'], categorical_features=cat_features,
verbose=True)
exp = explainer.explain_instance(arTest[10], model1.predict, num_features=6)
exp.show_in_notebook()
Opis zmiennych wybranych przez explainer i możliwe wartości:
- PREGNT31 - (-1,1,2) czy w okresie referencyjnym jest/była w ciąży: -1 - nie dotyczy, 1 - tak, 2 - nie.
- PCS42 - sumaryczna ocena zdrowia fizycznego (-1, wartość) - zmienna ciągła, jeśli wartość = -1 tzn., że ten wskaźnik nie mógł zostać obliczony: im wyżej tym lepiej.
- DIABDX - cukrzyca (-1,1,2)
- ACTLIM31 - ograniczenia związane z pracą/pracami domowymi (-1,1,2)
- CANCERDX - nowotwór (-1,1,2)
- RTHLTH31 - samoocena stanu zdrowia (-1, 1-5): im wyższa tym gorzej
Uwagi:
Model jako najistotniejszą cechę wybrał fakt bycia przez pacjenta w ciąży (w okresie, którego dotyczy zebrany wywiad). Wydaje się to mieć sens, jeśli kobieta była w ciąży to na pewno korzystała z usług medycznych - cały okres okołociążowy obiftuje w wizyty i badania lekarskie. Inna sprawa, że ta wybrana obserwacja dotyczy mężczyzny.
Druga wskazana cecha to fakt, że wskaźnik określający ogólny stan zdrowia fizycznego nie mógł zostać wyliczony (lub jest bardzo niski?). Być może zmienna ta wymaga modyfikacji.
Następnie model wskazuje na fakt niezdiagnozowania u pacjenta cukrzycy => wydatki na OM mniejsze.
Wyjaśnienie mówi też, że brak ograniczeń w wykonywaniu prac zarobkowych/domowych ujemnie wpływa na wydatki na opiekę medyczną.
Fakt niezdiagnozowania nowotworu ujemnie wpływa na przewidywane wydatki na OM.
Trudno jest zinterpretować wskazania modelu w przypadku samooceny stanu zdrowia (RTHLTH31). Wyjaśnienie mówi, że jeśli wartość ta jest różna od 4, to wydatki są generalnie mniejsze.
df_raw['RTHLTH31'].value_counts()
dfTest.iloc[obs_num_1,:]
df_raw['PHQ242'].hist(bins=2*7);
plt.plot()
several_obs = [23,32,42]
dfTest.iloc[several_obs,:]
for obs_id in several_obs:
explainer.explain_instance(arTest[obs_id], model1.predict, num_features=6).show_in_notebook()
Zmienne, które nie pojawiły się wcześniej:
- MCS42 - sumaryczna ocena zdrowia psychicznego (-1, wartość) - zmienna ciągła, im wyżej tym lepiej.
- WLKLIM31 - problemy z poruszaniem się (-1,1,2)
- EDRECODE- wykształcenie (-1,1,2,13,14,15,16) - im wyżej tym lepsze.
Widzimy, że wśród tych kilku wybranych obserwacji pojawia się zestaw powtarzających się cech i ich wpływ na model (jakościowo) jest taki sam. Najistotniejsza pozostaje zmienna dotycząca ciąży. Wszystkie zmienne, które odpowiadają za fakt zdiagnozowania/niezdiagnozowania konkretnych chorób (kończące się na DX) wydają się sensownie wpływać na wyjścia modelu (jeśli ktoś ma/miał zdiagnozowaną daną chorobę to wydatki rosną). Dodatkowo pojawiają się zmienne związane z samopoczuciem i oceną stanu zdrowia, których interpretacja znów jest utrudniona.
Jeden z modeli wskazał, że fakt iż osoba na pewno nie jest doktorem/magistrem (EDRECODE_16=0) powoduje, że wydatki na OM są mniejsze.
Uwaga ogólna - kodowanie jedynkowe utrudnia odczyt takich wykresów.
model2 = models['ANN']
model2
ar = np.asarray(dfTest)
model2.predict(arobs)
exp1 = explainer.explain_instance(arTest[10], model1.predict, num_features=6)
exp2 = explainer.explain_instance(arTest[10], model2.predict, num_features=6)
exp1.show_in_notebook()
exp2.show_in_notebook()
Zmienne, które nie pojawiły się wcześniej:
- K6SUM52 - ogólna ocena odczuć (mentalnych) (-1, 0-24): im wyżej, tym gorzej.
- Age31X - wiek
- RACE3 - rasa - (0,1): 0 - white, 1 - black
df_raw['RACE3'].value_counts()
for obs_id in several_obs:
explainer.explain_instance(arTest[obs_id], model2.predict, num_features=6).show_in_notebook()
Porównanie modeli dla obserwacji nr 10 i kilku losowych:
- Model ANN wskazał zupełnie inne cechy od modelu XGB
- Model ANN wydaje się mniej stabilny - dla różnych obserwacji wskazuje zupełnie inne cechy. Widać więcej parametrów ciągłów, podczas gdy model XGB preferował parametry kategoryczne.